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Abstract 

Background: Root length and its architecture govern the adaptability of plants to various stress conditions, including 
drought stress. Genetic variations in root growth, length, and architecture are genotypes dependent. In this study, we 
compared the drought-induced transcriptome of four genotypes of Gossypium herbaceum that differed in their drought 
tolerance adaptability. Three different methodologies, namely, microarray, pyrosequencing, and qRT-PCR, were used for 
transcriptome analysis and validation. 

Results: The variations in root length and growth were found among four genotypes of G.herbaceum when exposed to 
mannitol-induced osmotic stress. Under osmotic stress, the drought tolerant genotypes Vagad and GujCot-21 showed a 
longer root length than did by drought sensitive RAHS-14 and RAHS-IPS-187. Further, the gene expression patterns in the 
root tissue of all genotypes were analyzed. We obtained a total of 794 differentially expressed genes by microarray and 
104928 high-quality reads representing 53195 unigenes from the root transcriptome. The Vagad and GujCot-21 respond to 
water stress by inducing various genes and pathways such as response to stresses, response to water deprivation, and 
flavonoid pathways. Some key regulatory genes involved in abiotic stress such as AP2 EREBP, MYB, WRKY, ERF, ERD9, and 
LEA were highly expressed in Vagad and GujCot-21. The genes RHD3, NAP1, LBD, and transcription factor WRKY75, known 
for root development under various stress conditions, were expressed specifically in Vagad and GujCot-21. The genes 
related to peroxidases, transporters, cell wall-modifying enzymes, and compatible solutes (amino acids, amino sugars, 
betaine, sugars, or sugar alcohols) were also highly expressed in Vagad and Gujcot-21. 

Conclusion: Our analysis highlights changes in the expression pattern of genes and depicts a small but highly specific set 
of drought responsive genes induced in response to drought stress. Some of these genes were very likely to be involved 
in drought stress signaling and adaptation, such as transmembrane nitrate transporter, alcohol dehydrogenase, pyruvate 
decarboxylase, sucrose synthase, and LEA. These results might serve as the basis for an in-depth genomics study of 
Gossypium herbaceum, including a comparative transcriptome analysis and the selection of genes for root traits and 
drought tolerance. 
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Background 

During this century, water limitations are likely to in- 
crease in agriculture, and hence, the development of 
crop varieties with better water use efficiency (WUE) is 
of prime importance in agriculture. The root signifi- 
cantly contributes to the management of water stress 
and the adaptability of plants to stress [1]. Thus, 
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recently, a significant number of studies focusing on 
genetic variation in root development and its response 
to drought and other abiotic stress were reported in 
wheat [1], soybean [2], and maize [3,4]. Root length 
depicts the adaptability for drought stress in Oryza 
sativa [5,6]. Larger root systems provide the ability for 
plants to compete for nutrients, give them support for 
water uptake, and enable them to survive in water defi- 
ciency [7,8]. Root system development and its 
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architectures is determined by intrinsic genetic proper- 
ties; modulated by various environmental factors, nutri- 
ents, and water availability; and presumably governed by 
many genes [9]. Several genes that control root architec- 
ture and development have been identified in Arabidop- 
sis thaliana, Zea mays, and Oryza sativa in recent 
genetic studies [10-12]. The list of genes includes 
LBD16/ASL18, LBD29/ASL16, and CRL1/ARL1, which 
are reported to regulate lateral root formation and help 
in root gravitropism [10]. Further, RHD3 was reported 
to be responsible for the formation and elongation of 
root hairs [13]. The transcription factors HD-ZIP and 
PLETHORA (PLT) were primarily responsible for the 
molecular origin of the root and cell-type differentiation 
[14,15]. Similarly, other transcription factors, namely 
WRKY75, ZAT6 [16-18], and bHLH32 [19], have been 
reported to play a critical role in root development and 
growth under phosphate deficiency. The root plays a 
vital role in maintaining the physiological condition of 
plants during water stress. The transpiration rate and 
stomatal conductance of plants are reduced during 
drought stress, and they are stimulated by chemical and 
hormonal signaling before hydraulic signaling in the 
roots [20]. Various signaling molecules such as auxin 
and cytokinin are produced in the roots and play a cru- 
cial role in shoots during the drought stress in plants 
[21]. The auxin helps in cell division, tropisms, vascular 
differentiation, and maintenance of the root meristem 
[22]. Simultaneously, in response to external stimuli 
such as drought, it changes the transcriptional programs 
of specific cells types and transmits the spatial informa- 
tion in different plant organs [23]. Cytokinin signaling 
is mediated through histidine kinase pathways in root 
tissues under the drought condition [24]. Cytokinin reg- 
ulates the pattern formation and differentiation of the 
vascular tissue of roots during root development [25]. 
Various signal transduction cascades mediating hormo- 
nal signaling in response to nutrient deprivation and 
drought stress and their transcriptional responses have 
been elucidated from microarray and transcriptome 
sequencing studies [26]. Drought stress creates a redox 
imbalance in plant cells, resulting in a cascade of 
changes in various signaling pathways, cell wall- 
membrane integrity, and produces various antioxidant 
enzymes [27]. Many of these enzymes, including super- 
oxide dismutases (SODs), ascorbate peroxidases (APXs), 
catalases (CATs), glutathione peroxidases, and peroxire- 
doxins, help in the maintenance of the redox balance of 
cells under the stress condition [28]. The roles of various 
membrane transporters have been characterized under 
drought and salt stress, which provide osmotic balance 
and ionic homeostasis in plant cells [29]. 

The diploid (2n = 26) species of cotton, namely G. her- 
baceum and G.arboreum, are indigenous in Asia and 



Africa and are commonly referred to as desi cottons in 
India [30,31]. The Gossypium herbaceum (Al-genome) 
and Gossypium arboreum (A2-genome) have an inherent 
ability to sustain drought stress [32]. In drought-prone 
areas in Asia, these two species are preferentially culti- 
vated due to their ability to withstand drought [32]. We 
recently published a detailed physiological investigation 
and conducted genome-wide expression studies on 
leaves of drought-tolerant Vagad and drought-sensitive 
RAHS-14 genotypes [31]. We identified the various 
genes and inherent biochemical pathways that interact 
in drought-tolerant genotypes to provide them with an 
adaptive advantage in the stress condition [31]. Here, we 
report the expression profiling of the root of four geno- 
types of Gossypium herbaceum, namely Vagad, GujCot- 
21 (drought tolerant), RAHS-14, and RAHS-IPS-187 
(drought sensitive) using Affymetrix cotton expression 
array and GS-FLX pyrosequencing in response to 
drought stress. We identified many stress-responsive 
genes and biochemical pathways that were differentially 
expressed in roots during the drought condition in the 
selected genotypes. Further, we could assign these 
root-specific genes to different root zones based on 
the Genevestigator tool. These differentially expressed, 
tissue-specific genes in the contrasting genotypes re- 
flect the probable mechanistic details of the adaptabi- 
lity of tolerant genotypes to water stress. 

Methods 

Measurement of root length 

Cotton (G.herbaceum) seeds of Vagad, GujCot-21, 
RAHS-IPS-187, and RAHS-14 were sterilized, aseptically 
immersed in water for a day at 30°C, and then placed for 
germination in a moist petri dish under the following 
condition: 28°C/25°C as day and night temperature, 12 h 
of light and dark period alternatively, and a relative hu- 
midity of 80%. After 3 days, the properly germinated 
seeds were transferred to paper rafts, which were sub- 
merged (3/4 th ) in Hoaglands media containing a differ- 
ent percentage of mannitol (4%, 6%, and 8%), and only 
Hoaglands media was considered as a control [31]. The 
seedlings were then allowed to grow until they displayed 
healthy growth in the control condition; next, the root 
lengths of the grown seedlings of all genotypes were 
photographed and measured. 

Drought treatments and growth condition 

Cotton (G.herbaceum) seeds of four genotypes, namely 
Vagad, RAHS-14, GujCot-21, and RAHS-IPS-187, were 
germinated with six seedlings each in an earthen pot 
in the soil (one plant/pot) and grown in a greenhouse 
(16-h-light/8-h-dark cycles) at 28°C to 30°C with a rela- 
tive humidity of 50-60% and a photosynthetically active 
radiation of 900 [imol m- 2 s- 1 . After 4 weeks, plants 
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containing six leaves were selected for the experiments 
and drought stress was imposed by withholding water 
for seven days (soil moisture below 30%), and control 
plants were watered regularly [31]. After 7 days of 
drought treatment, drooping effects on plant leaves be- 
came prominent; the plants were uprooted, and the 
roots of both control and drought treated plants were 
homogenized in liquid nitrogen before proceeding with 
RNA isolation. 

Isolation of RNA and real-time PCR 

Total RNA was isolated in three independent biological 
replicates of each genotype using the Spectrum Plant 
Total RNA Kit (Sigma- Aldrich) according to the manu- 
facturers instructions and eluted with nuclease-free 
water. After DNasel (Ambion) treatment, they were 
quantified, and checked for integrity using a Bioanalyzer 
(Agilent, Inc., Palo Alto, CA, USA). For the analysis of 
target genes by RT-PCR, cDNA was synthesized from 
the total RNA of Vagad and RAHS-14 by reverse tran- 
scriptase using the cDNA synthesis kit (Invitrogen) 
according to the manufacturers instructions. The RT- 
PCR was performed in replicates with 1 \A cDNA using 
the SYBR Green Master Mix (Applied Biosystems) with 
an ABI 7500 sequence detection system as prescribed in 
the manufacturers protocol (Applied Biosystems). The 
gene-specific and ubiquitin (internal control) primers 
were designed using PRIMER3 (http://frodo.wi.mit.edu/ 
primer3/input.htm). The primer details are listed in 
Additional file 1. The relative quantification method 
(AACT) was used for the quantitative evaluation of the 
gene expression. 

Affymetrix gene chip experiment and microarray 
data analysis 

Three sets of biological replicates of Vagad and RAHS- 
14 under control and drought conditions (3x2x2) were 
taken independently, constituting a total of twelve cot- 
ton chips for analysis. For Affymetrix gene chip analysis 
(Affymetrix, Santa Clara, CA, USA), 1 \ig total RNA of 
root tissue was used for making biotin-labeled cRNA 
targets, hybridization. Synthesis of cDNA, cRNA, and its 
fragmentation, hybridization, washing, staining, and 
scanning were conducted according to the gene chip 
standard protocol (Eukaryotic Target Preparation; Affy- 
metrix). The signal intensity of each probe set on the 
cotton gene chip was analyzed using Affymetrix GCOS 
software, and the target mean value (TGT) was scaled as 
being 500 for each chip. Students t-test analysis and 
log 2 -transformed signal ratio of each probe set were car- 
ried out by the Array Assist Software 5.2.2 (Agilent 
Technologies, Santa Clara, CA, USA). Differentially 
expressed genes with a detection p value < 0.05 were 
considered present in three biological replicate 



experiments. Gene expression data analyses were com- 
pleted using a filtered RMA expression value. The anno- 
tation of the probe set in the Affymetrix cotton genome 
array was mapped with the locus ID of Arabidopsis 
TAIR10 version by BLAST. Based on the annotation, the 
expressed genes were analyzed. Gene ontology analysis 
was performed for the functional categorization of dif- 
ferentially expressed genes using agriGO tool (http:// 
bioinfo.cau.edu.cn/agriGO/), and the p-values were cor- 
rected by applying the false discovery rate correction to 
control falsely rejected hypothesis during the analysis. 
The FDR corrected p-value of <0.05 was assumed as the 
cutoff value. Microarray gene expression data used in 
this study were MIAME compliant and deposited in 
NCBI-GEO with accession number GSE36249. 

Double strand cDNA library construction and GS-FLX 
pyrosequencing 

Total RNA (3 \ig) from the root tissue of GujCot-21 and 
RAHS-IPS-187 was reverse transcribed using a T7-OHgo 
(dT) Promoter Primer in the first-strand cDNA synthesis 
(Affymetrix). After RNase H-mediated second-strand 
cDNA synthesis, the double-stranded cDNA was 
enriched and served as a template in the subsequent 
in vitro transcription (IVT) reaction (Affymetrix). The 
IVT reaction was carried out in the presence of T7 RNA 
Polymerase (Affymetrix). The cRNA (3 \ig) was reverse 
transcribed in the first-strand cDNA synthesis step by 
using a random hexamer primer, followed by RNase 
H-mediated second-strand cDNA synthesis in replicates. 
The replicate samples were pooled and purified by a 
column (QIAquick PCR purification kit) and further 
used for GS-FLX pyrosequencing. 

Emulsion based clonal amplification 

Double-strand cDNAs obtained from Gujcot-21 and 
RAHS-IPS-187 were used for GS FLX library prepar- 
ation. Approximately 5 \xg of double-strand cDNA was 
sheared by nebulization at 206 kPa for 2-4 min. The 
fragmented cDNA were amplified in aqueous droplets 
that were made through the creation of a PCR reaction 
mixture in emulsion oil (Roche Diagnostics). The dro- 
plets acted as separate microreactors in which parallel 
DNA amplifications were performed, yielding approxi- 
mately 10 7 copies of a template per bead. After PCR, the 
emulsion was broken to release the beads containing the 
amplified DNA template. The beads carrying the tem- 
plates were enriched and deposited by centrifugation 
into the open wells of a 70 x 70 mm 2 optical picotiter 
plate for sequencing. 

Assembly and annotation of transcriptome 

All sequence analyses were conducted using a publicly 
available software, R package (www.R.proj ect.org). The 
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pyrosequencing reads were assembled after quality con- 
trol using software version 2.5 newbler, and the assembly 
was performed using a 40 bp overlap length and a 90% 
identity. Transcripts annotation was performed by 
BLASTX using the Basic Local Alignment Tool against 
the NR database from NCBI (http://www.ncbi.nlm.nih. 
gov/); from TAIR (http://www.arabidopsis.org/) and the 
BLASTN for EST cotton database available at NCBI. 
The GS-FLX sequence reads discussed in this article can 
be found in the Genebank (http://www.ncbi.nlm.nih.gov/ 
genbank) of the National Center for Biotechnology with 
accession number SRA029162.1. 

Digital expression analysis 

For the digital expression analysis, the reads for both li- 
braries (GujCot-21 and RAHS-IPS-187) were tagged and 
pooled to form a large dataset of 104928 reads. These 
reads were assembled using the newbler assembler at an 
overlapping of 40 bp and a 90% identity. These reads 
were assembled into 2080 contigs. We calculated TPM 
value and R value using the R statistics for supercontigs 
and and those with R value > 3 and Fold change > 2 were 
considered significantly differentially expressed contigs 
[33]. These filtered contigs were annotated using 
BLASTX against the NCBI-NR database, the TAIR, and 
BLASTN was used for the cotton EST public database. 
The Genevestigator (https://www.genevestigator.com) 
software was used for the tissue-specific expression ana- 
lysis of differentially expressed genes. 

GO and KOBAS analysis 

The functional classification of unigenes was performed 
by assigning gene ontology annotation codes. Only those 
unigenes that expressed a significant Blast result against 
the TAIR database were used for GO annotation. The 
GO annotation for level three was extracted for each 
library and used for further analysis. The biochemical 
pathway assignments were carried out according to 
kobas (http://kobas.cbi.pku.edu.cn/home.do). 

Results 

Root growth determined drought tolerance properties in 
cotton 

Different genotypes of G.herbaceum were screened for 
apparent differences in the root growth under control 
and osmotic stress conditions. The four genotypes, 
namely Vagad, GujCot-21, RAHS-14, and RAHS-IPS- 
187, exhibited contrasting differences in their root 
structure under control and osmotic stress conditions 
(Figure 1). Vagad and GujCot-21 have a longer root 
length as compared with RAHS-14 and RAHS-IPS-187 
in the control condition. The difference was most pro- 
nounced in osmotic stress treatment by different con- 
centrations of mannitol, which showed a significantly 



longer root length in Vagad and GujCot-21 at 6% of 
mannitol concentration as compared with a root length 
of RAHS-14 and RAHS-IPS-187. Further, at 8% of man- 
nitol concentration, Vagad and GujCot-21 showed 
stunted growth of the root, but in RAHS-14 and 
RAHS-IPS-187, the development of the root was com- 
pletely abolished (Figure 1). For further molecular 
analysis of the root of these genotypes, Vagad and 
RAHS-14 were analyzed by microarray, and GujCot-21 
and RAHS-IPS-187 were analyzed by pyrosequencing. 

Comparative gene expression analysis of root tissue of 
Vagad and RAHS-14 by Affymetrix cotton gene chip 

Expression profiling experiments conducted on cotton 
root samples from drought-tolerant (Vagad) and sensi- 
tive (RAHS-14) genotypes were compared by a commer- 
cial Affymetrix cotton gene chip, which is represented 
by 21,854 transcripts of cotton. The differentially 
expressed transcripts of the root from Vagad and 
RAHS-14 under drought stress and control conditions 
were determined. The genes that expressed equal or 
greater than two fold at a p-value < 0.05 were assigned 
as differentially expressing genes for further analysis. 

Differentially expressed (DE) gene analysis 

Our analysis identified 156 and 538 transcripts as being 
differentially up-regulated at the FC > 2, p < 0.05 signifi- 
cance level in RAHS-14 under water and drought condi- 
tions, respectively (Additional file 2 and Additional file 
3). Similarly, 165 and 256 transcripts were identified as 
being differentially up-regulated at the FC > 2, p < 0.05 
significance level in Vagad under watered and drought 
conditions, respectively (Additional file 4 and Additional 
file 5). Out of 538 differentially expressed genes in 
RAHS-14 under the drought condition, 80 genes exhib- 
ited an expression more than five fold as compared with 
Vagad root in the drought condition. These genes in- 
clude cyclopropane fatty acid synthase, transcription 
regulator NOT family transcription factor, curculin-like 
lectin proteins, pectin esterase, metal ion transmem- 
brane transporter, cytochrome p450, brassinosteroid 
insensitive 1- associated receptor kinase, and heat-shock 
protein family (Additional file 3). Many of the genes in 
RAHS-14 (probe sets 538) represent stress-responsive 
genes as per their GO annotation and were previously 
shown to be involved in abiotic stress response. Simi- 
larly, 156 genes that expressed more than two fold in 
RAHS-14 under the water condition include senescence- 
associated protein (FC > 50), zinc finger protein, 
dehydration-associated proteins, AP2 domains contain- 
ing transcription factor, and NAC domain transcription 
factor observed a higher expression in RAHS-14 
(Additional file 2 ). The 165 genes uniquely up-regulated 
in Vagad under the water condition include metal 
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Figure 1 Effect of osmotic stress induced by different percentages of mannitol in the hydroponics condition on the development of 
the root of Vagad, GujCot21, RAHS-14, and RAHS-IPS187. 



binding proteins (FC > 200), osmotin-like protein, cyto- 
chrome, and various other genes and exhibited more 
than two fold expression (Additional file 4). While in 
drought stress, the 256 up-regulated genes include Xylo- 
glucan endotransglycosylase, serine proteases, lipid- 
binding protein, glycine-rich protein, sodium/potassium 
proton exchanger, and aquaporins proteins (Additional 
file 5). 

Gene ontology enrichment and Kyoto Encyclopedia of 
Genes and Genomes (KEGG) pathway analysis of 
differentially expressed genes 

To identify the possible biological pathways that govern 
the responses of differentially expressed genes, gene 
ontology (GO) and KEGG pathway analysis was carried 
out for Vagad and RAHS-14 under water and drought 
stress. The GO annotation mapping of RAHS-14 



revealed that the genes assigned to GO terms for the 
thiamine metabolic process (6 genes) and response to 
auxin (5 genes) were exclusively enriched in the water 
condition (Figure 2). However, under drought stress, the 
genes cover a broad range of GO categories that repre- 
sents response to abiotic stimulus, response to biotic 
stimulus, sucrose metabolic process, response to heat, 
cell death, and peptidase activity were enriched only 
under the drought condition in RAHS-14 (Figure 2). 
Various other stress-responsive GO terms significantly 
present under the drought condition were transferase ac- 
tivity, cellular ketone metabolic process, oxidoreductase 
activity, response to stress, transition metal ion binding, 
amine metabolic process, kinase activity, ATP binding, 
macromolecular complex, and integral to membrane 
(Figure 2). A small set of GO terms were over- 
represented under the water condition, such as ligase 
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activity, protein binding, transcription factor activity, and 
chloroplast related genes in RAHS-14 (Figure 2). Other 
GO terms with a nearly similar expression level under 
both water and drought conditions in RAHS-14 were 
response to hormone stimulus, transcription regulator 
activity, and DNA binding activity in RAHS-14 
(Figure 2). 

Similarly, GO annotation mapping in Vagad under the 
drought condition showed that various differentially 
expressed genes assigned to GO terms were involved in 
nucleic acid binding, response to abiotic stimulus, cellu- 
lar carbohydrate metabolic process, and ATP binding 
(Figure 3). Likewise, the GO terms covered a broad range 
of stress and membrane integrity-related categories that 
were significantly enriched under drought as compared 



with the water, which includes the response to stress, 
hydrolase activity, macromolecular complex, regulation of 
the biological process, integral to membrane, intrinsic to 
membrane, transporter activity, cellular metabolic process, 
and so on (Figure 3). A single GO term response to biotic 
stimulus was present differentially under the water condi- 
tion (Figure 3). Under the water condition, some interest- 
ing GO terms were over-represented in Vagad, which 
include oxidoreductase activity, ion binding, metal ion 
binding, catabolic process, cellular amine metabolic 
process, cellular nitrogen compound metabolic process, and 
so on (Figure 3). When subjected to drought stress, as an 
adaptative mechanism, Vagad underwent re-programming 
of a large number of genes; hence, almost all GO terms 
exhibited vast differences at the expression level 
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Figure 2 GO annotation of differentially expressed up-regulated genes in RAHS-14 during control and drought conditions. 
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Figure 3 GO annotation of differentially expressed up-regulated genes in Vagad during drought and control conditions. 



Microarray validation: quantitative RT-PCR analysis 

To verify the microarray results, quantitative PCR assays 
were carried out for the 12 selected genes. Three 
genes, namely Ghi.6528.1.Sl_at, GraAffx.8742.1.Sl_at, and 
GraAffx.33038.1.Sl_s_at, were identified as being up- 
regulated in response to drought stress, and three genes, 
namely squalene monooxygenase (Ghi.482.1.Al_s_at), 
GhiAffx.46297.1.Sl_s_at, and GhiAffx.43814.1.Sl_at, were 
identified as being down-regulated in Vagad in microarray 
data analysis and were selected for further validation 
(Figure 4A). Similarly, three genes, namely GhiAffx.60321.1. 
Sl_x_at, Ghi.6435.1.Al_at, and GhiAffx.31372.1.Sl_at, 
were identified as being up-regulated in response to 
drought stress, and three genes, namely senescence- 
associated protein (Ghi.743.1.Sl_at), Ghi.l092.3.Sl_s_at, 
and Ghi.3284.1.Sl_s_at, were identified as being down- 
regulated in RAHS-14 and were selected for further valid- 
ation by using quantitative RT-PCR (Figure 4B). The 
quantitative RT-PCR analysis demonstrated that the 
up-regulated genes in Vagad showed 5- to 20-fold higher 



expressions in response to drought stress and that the 
down-regulated genes showed 5- to 30-fold repression in 
response to drought. The results were similar in case of 
the genes selected for RAHS-14. Thus, quantitative RT- 
PCR results agree with microarray analysis and, hence, 
validate the microarray data. 

Gene expression analysis of tolerant (GujCot-21) and 
sensitive genotype (RAHS-IPS-187) under drought stress 
by pyrosequencing of transcriptomes 

Our initial analysis suggested that GujCot-21 was yet 
another drought-tolerant genotype, and RAHS-IPS-187 
was a drought-sensitive genotype in our studies (Fig- 
ure 1). Thus, the genotypes GujCot-21 (drought tolerant) 
and RAHS-IPS-187 (drought sensitive) of G. herba- 
ceum were taken for further analysis by transcriptome se- 
quencing of their root tissue to investigate gene 
expression patterns under drought stress. The 55620 and 
49308 sequencing reads were obtained from GujCot-21 
and RAHS-IPS-187, respectively (Table la). The average 
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Figure 4 Validation of microarray differentially expressed genes by quantitative PCR. (A) Vagad (B) RAHS-14. 



lengths of the assembled contigs and singletons were 
nearly 481 bp and 237 bp, respectively (Table la). The 
average depth of the contigs in GujCot-21 and RAHS- 
IPS -187 was about five reads per contig. For differential 
expression analysis of the genes in GujCot-21 and RAHS- 
IPS -187, the reads from both genotypes were tagged and 
pooled to form one large dataset that was assembled into 
contigs by using Roches GS-assembler (overlap length 
>40 bp and 90% identity). The 104928 reads were 
clustered into 2664 contigs and 50,531 singletons with an 
average length of 508 bp and 231 bp, respectively, 
for the expression analysis (Table la). The significant 
changes in gene expression as transcript per million 



(TPM) were calculated using the R statistics (R value > 3) 
(Additional file 6). 

Annotation of root transcriptome of Gossypium 
herbaceum 

The resulting contigs and singletons were queried 
against the NBCI NR, TAIR database for annotation 
using the BLASTX program at a stringency of e-value of 
10E-5 and greater than 50% overlap. Each library was 
annotated with these databases separately. Out of the 
total contigs and singletons in GujCot-21, 1011 contigs 
and 17081 singletons were annotated with the TAIR 
database (Table lb, Additional file 7), 1056 contigs and 



Table 1 Summary of transcriptome sequencing data 



(a) 




Reads 


Bases 


Contigs 


Singleton 


Av. Contig length 


Av. S. length 


GujCot-21 


55, 620 


13, 020, 140 


1,281 


30, 501 


481.7 bp 


237.8 bp 


RAHS-IPS 187 


49, 308 


11, 199, 207 


858 


30, 776 


532.9 bp 


228.6 bp 


Supercontigs 


1, 04, 928 


24,219, 347 


2, 664 


50, 531 


508.7 bp 


231.7 bp 


(b) 




TAIR 




NR 




Cotton EST 






Contigs 


Singleton 


Contigs 


Singleton 


Contigs 


Singleton 


GujCot-21 


1, 011 


17, 081 


1, 056 


18, 392 


779 


14,476 


RAHS-IPS- 187 


682 


16, 856 


713 


9, 885 


488 


12, 383 


Super-Contigs 


2, 080 


26, 838 


2, 163 


17204 


1, 537 


20, 552 


(c) 


Contigs 


TAIR 


NR 


Cotton_EST 


Total hits 


No- hits 




GujCot-21 


1, 011 


1, 056 


779 


1266 


15 




RAHS-IPS-187 


682 


713 


488 


855 


2 




Supercontigs 


2, 080 


2, 163 


1, 537 


2432 


232 





(a) Sequencing and assembly output (b) Annotation of contigs against NCBI-NR,TAIR and cotton EST (c) total number of annotated and novel contigs (no hit). 
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18392 singletons with the NR database, respectively 
(Additional file 8). In RAHS-IPS 187, 682 contigs and 
16856 singletons were annotated with TAIR, 713 contigs 
and 9885 singletons with the NR database, respectively 
(Table lb). For the digital expression analysis, the 
assembled data of GujCot-21 and RAHS-IPS-187 were 
also annotated with these databases. Approximately 
81.25% of the assembled data produced significant hits 
with the TAIR and NR database. A total of 2080 super- 
contigs and 26838 singletons were annotated with the 
TAIR, 2163 supercontigs and 17204 singletons with the 
NR, respectively. To find the common sequences be- 
tween already reported cotton ESTs and our unigenes, 
we queried the our dataset against all publicly available 
cotton ESTs, at criteria of e value of 10" 5 , and at least 
99% alignment of either the query or the subject. At 
these criteria, 779 and 488 contigs of GujCot-21 and 
RAHS-IPS-187, respectively, matched with cotton ESTs 
(Additional file 9). Only 15 sequences from GujCot-21 
and 2 fromRAHS-IPS-187 did not have any match with 
NR, TAIR and cotton EST database and hence were 
considered unique genes (Table lc). The number of 
contigs was assigned as unique and common in NCBI; 
the TAIR and cotton EST databases and presented in 
Additional file 10. A total of 15, 2 , and 232 contigs were 
uniquely present in all the 3 database of tolerant geno- 
types, sensitive genotypes, and supercontigs, respectively. 

Transcriptome analysis of cotton under drought stress 
reveals large number of drought induced genes and 
novel transcripts 

Differentially expressed genes between GujCot-21 and 
RAHS-IPS187 were calculated using the R statistics (R 
value > 3), and 2026 genes were differentially expressed 
in both the genotypes.For each of the contigs, the counts 
were converted to transcripts per million (TPM), which 
were transformed [log2 (fold change values)], and their 
differences were calculated for the fold change between 
GujCot-21 and RAHS-IPS-187 (Additional file 6). The 
quantitative profiling of the transcriptome reveals that 
1503 and 1160 genes were differentially up-regulated in 
GujCot-21 and RAHS-IPS-187 genotypes, respectively, 
as compared with each other. Interestingly, 135 (17.95%) 
and 96 (18.17%) genes were up-regulated in GujCot-21 
and RAHS-IPS-187 genotypes, respectively, and showed 
no hit in any database. Further, to obtain a global view 
of the cotton transcriptome and gene activity in two 
contrasting genotypes GujCot-21 and RAHS-IPS-187, 
differentially expressed genes were compared on the 
basis of their functional annotation. We found that the 
over-representation of the transcripts related to oxidore- 
ductases, ribosomal proteins, membrane transporter, cal- 
cium ion binding proteins, ATPases, dehydrogenases, 
heat shock protein, and various enzymatic process 



related to other activities in GujCot21, where as in case 
of RAHS-IPS-187 transcripts related to hydrolases, 
nuclear proteins and structural proteins like structu- 
ral constituent of cytoskeleton were over-represented 
(Additional file 6). There are several possible implica- 
tions for the enrichment of these functional classes of 
genes in specific genotypes. The relative abundance of 
oxidoreductases, membrane transporter, cell wall-related 
proteins, and heat shock protein in GujCot-21 may be 
directly related to their drought-tolerant behavior. 

Functional annotation of differentially expressed genes 
by GO and KOBAS 

In order to assess the differentially expressed genes with 
their associated biological pathways and to study the 
plant adaptations to limited water condition in the 
drought-contrasting genotype of cotton, several unique 
and common pathways affected by drought stress in 
GujCot-21 (1504 genes involved) and RAHS-IPS-187 
(1161 genes involved) were identified. The significant 
GO terms (p-value < 0.05) were characterized into bio- 
logical processes and molecular functions by agriGO 
(http://bioinfo.cau.edu.cn/agriGO/). The most significant 
categories, namely response to water deprivation (24%), 
response to jasmonic acid stimulus (21%), ethylene- 
mediated signaling pathway (10%), hyperosmotic re- 
sponse (7%), hyperosmotic salinity response (6%), water 
transmembrane transporter activity (6%), and so on, 
were the major abiotic stress-related GO categories ex- 
clusively present in GujCot-21 and, hence, gained atten- 
tion. Other GO terms exclusively present in GujCot-21 
include response to inorganic substance (30%), methio- 
nine metabolic process (9%), anatomical structure ar- 
rangement (8%), response to brassinosteroid stimulus 
(6%), negative regulation of signal transduction (5%), 
and so on. Also noteworthy were the categories of re- 
sponse to heat (22%), water-soluble vitamin biosynthetic 
process (8%), response to sucrose stimulus (8%), starch 
biosynthetic process (5%), regulation of stomatal move- 
ment (5%), and so on, which were exclusively present in 
RAHS-IPS-187 with interesting GO terms. The genes 
associated with 4 common metabolic pathways include 
response to reactive oxygen species (12% and 9% of 
associated genes), response to hydrogen peroxide (8% 
and 6% of genes), response to high light intensity (7% 
and 9% of genes), and lipid localization (5% and 7% of 
genes) in GujCot-21 and RAHS-IPS-187, respectively 
(Figure 5 and additional file 6). In addition, the KOBAS 
pathway analysis was also performed to determine the 
various specific pathways involved in the differential up- 
regulation of genes. Several biochemical pathways were 
mapped from gene sets of GujCot-21 and RAHS-IPS- 
187; among them, significant KOBAS -terms were fil- 
tered with a stringent criterion of a p-value < 0.05. 
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Except one pathway, all the other KOBAS terms fil- 
tered were exclusively present in either of the geno- 
types. Significantly, nine KOBAS pathways were 
specifically modulated in GujCot-21 and six in RAHS- 
IPS-187 (Additional file 6). A large number of genes 
from GujCot-21 (91 genes) were mapped to the ribo- 
some; ath03010, making highest a p-value of 6.10E-03, 
and 33 genes form RAHS-IPS-187 with a p-value of 
9.52E-05. Spliceosome; ath03040 is the second-most 
significant in GujCot-21 and includes 17 associated 
genes (p-value - 0.006). Several other interesting 
and important pathways such as Glycolysis / Gluco- 
neogenesis (athOOOlO; P-value - 0.024017), Flavonoid 
biosynthesis (ath00941; P-value - 0.038), Fatty acid me- 
tabolism (ath00071; P-value - 0.043), Circadian rhythm 
- plant (ath04712; P-value - 0.044), Terpenoid back- 
bone biosynthesis (ath00900, P-value - 0.045), Histi- 
dine metabolism (ath00340; P-value - 0.052), and 
beta-alanine metabolism (ath00410; P-value - 0.054) 
were involved and synthesized in response to drought 
stress in GujCot-21 (Figure 6A). A pathway analysis of 
the genes present in RAHS-IPS-187 revealed plant 



hormone signal transduction (ath04075; P-value - 5.8 IE- 
05) to be the most significantly represented. It was fol- 
lowed by Phagosome (ath04145; P-value - 0.001), Endo- 
cytosis (ath04144; P-value - 0.005), Ubiquitin-mediated 
proteolysis (ath04120; P-value - 0.029), and Protein proces- 
sing in endoplasmic reticulum (ath04141; P-value - 0.033), 
which were biochemical pathways affected by drought 
stress in RAHS-IPS-187 (Figure 6B). These biochemical 
pathways present in GujCot-21 and RAHS-IPS-187 indi- 
cate the complexity of the metabolic changes involved in 
plant responses to drought stress. 

Higher expression of antioxidant enzymes may contribute 
to drought tolerance in cotton 

We identified the transcripts that encode various anti- 
oxidant enzymes which were predominantly expressed 
in GujCot-21, and, hence, these genes could be consid- 
ered probable components of the complex gene net- 
works involved in the drought condition. Hydrogen 
peroxide (H 2 0 2 ) and superoxide are the major reactive 
oxygen species (ROS) that are produced in plant cells 
under biotic or abiotic stresses [34]. Plants have 
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Figure 5 Gene ontology of differentially expressed contigs (fold change > 2) in drought-tolerant (GujCot-21) and sensitive genotypes 
(RAHS-IPS-187). MF-molecular function, BP-Biological process. 
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Figure 6 Functional enrichment of differentially expressed contigs (fold change > 2) by KOBAS (A) Tolerant genotype (GujCot-21) 
(B) Sensitive genotype (RAHS-IPS-187). 



developed several strategies to minimize the oxidative 
damage of ROS [28]. For instance, hydrogen peroxide is 
eliminated by catalases (CAT) and ascorbate peroxidases 
(APX), and superoxide is eliminated by superoxide dis- 
mutase (SOD). Cytosolic glutathione reductases (GR) 
also play a role in the detoxification of H 2 0 2 [35]. 
Microarray differentially expressed genes (FC > 2) of 



Vagad and RAHS-14 fetched 13 antioxidants genes, and 
GujCot-21 contigs from transcrip tomes contributed by 
way of 23 antioxidant-related transcripts. The specific 
expression of these 36 antioxidants was further classified 
in three categories, namely, superoxide dismutase, per- 
oxidase, and Glutathione S -transferase. Gujcot21 had 7 
contigs encoding superoxide dismutase, 12 contigs of 
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peroxidase, and 6 contigs related to Glutathione S -trans- 
ferase. The Cu/Zn superoxide dismutase activity showed 
a significant 6.75-fold higher expression in GujCot-21 
(contig 01875) (Table 2). Among several Glutathione 
S-transferases, the gene encoding GST group 8 was 
predominantly expressed in GujCot-21, and its activity 
was reported to play a major role in oxidative stress 
metabolism [36]. Out of 13 probe sets, 7 from RAHS- 
14 were distinguished into glutamate (6 probe sets) 
and peroxidase (1 preobe set), and the remaining 6 
probe sets from Vagad were grouped into peroxidase 
(4 probe sets) and glutathione transferase (2 probe sets). 
The probe set id GhiAffx.53252.1.Sl_s_at encoded glu- 
tamate synthase with highest 5.7-fold expression in 
RAHS-14. The 5.4- fold induction of glutathione trans- 
ferase (GhiAffx.36353.1.Sl_at) was highest in Vagad 
(Table 2). We found that most of the antioxidant 
enzymes were highly expressed by drought stress in 
GujCot-21, which implies that these antioxidant 
enzymes might have played important roles in plant 
tolerance to drought stress, whereas RAHS-IPS-187 
was susceptible to drought stress. 

Different regulations of transporter and cell wall-related 
genes between the drought-tolerant and 
sensitive genotypes 

Previous findings report that the membrane transpor- 
ters, which function in cellular transport processes for 
the maintenance and re-establishment of homeostasis in 
the plant cytoplasm, were induced with various forms of 
abiotic stress [37]. Both RAHS-14 and Vagad showed 
that a large number of genes belonging to the membrane 
transporter family significantly increased in response to 
drought stress (Table 3a). However, the most contrasting 
difference between Vagad and RAHS-14 was that in 
Vagad, many membrane transporters belonging to the 
ABC transporter family were expressed in response to 
drought stress in roots. Similarly, in case of RAHS-14, 
the transporters belonging to the ATPase super family 
were predominantly expressed in response to the 
drought stress (Table 3a). Further, a previous study 
revealed that the expression of the genes involved in cell 
wall metabolism was generally repressed during osmotic 
stress in Arabidopsis [38]. We notably observed that 22 
out of 28 genes belonging to cell wall biogenesis were 
repressed in both the genotypes; however, the expression 
of these genes was significantly higher in RAHS-14 
(Table 3b). For example, the expressions of two genes 
(Ghi.l0072.3.Al_at and GhiAffx.60321.1.Sl_x_at) that 
putatively encode glycine-rich protein and pectinesterase 
were largely unaffected during drought stress in Vagad, 
but they were significantly down-regulated in RAHS-14 
(Additional files 3 and Additional file 5). A comparative 
transcriptome analysis of drought-tolerant genotype 



(GujCot-21) and sensitive genotype (RAHS-IPS-187) 
also showed a similar expression pattern of the mem- 
brane transporter-related genes as had been obtained in 
a microarray analysis of Vagad (tolerant) and RAHS14 
(sensitive). Transcriptomic analysis showed approxi- 
mately 3.2% of differentially expressed genes in GujCot- 
21 encoding membrane transporter-related proteins. 
This result indicates that a comparatively less perturb- 
ation in the expression of the genes involved in cell wall 
biology will likely lead to better cell wall protection dur- 
ing drought stress in Vagad and GujCot-21 than in 
RAHS-14 and RAHS-IPS-187, thereby resulting in the 
better drought tolerance of Vagad and GujCot-21. Cell 
wall biogenesis plays an important role in root cell ex- 
pansion and, in turn, the root architecture; therefore, we 
next examined in our data the genes belonging to cell 
wall biogenesis. We found that the number of genes re- 
sponsible for cell wall biogenesis significantly increased 
in RAHS-14 and Vagad in response to drought stress 
(Table 3b). The genes representing the enzymes for cell 
wall biogenesis, namely cellulose synthase, glucosyl syn- 
thase, and epimerase, were found to be expressed in 
both Vagad and RAHS-14 in drought stress (Table 3b). 
However, in case of Vagad, we identified that xyloglucan: 
xyloglucan transferase, pectinesterase, lyases, and expan- 
sin genes were uniquely expressed in Vagad under 
drought stress as compared with RAHS-14. A higher 
expression of these genes further helps in cell wall bio- 
genesis and, hence, probably provides better drought 
tolerance in Vagad. 

Transcription factors in the drought stress response in 
roots of G.herbaceum 

To identify the putative transcription factors in the dif- 
ferentially expressed transcripts of all 4 genotypes, the 
cotton Affymetrix probe ID and pyrosequencing contigs 
were annotated with the corresponding homologous 
locus IDs of Arabidopsis and queried to the AGRIS data- 
base (http://arabidopsis.med.ohio-state.edu/AtTFDB/). A 
total of 138 genes encoding various related transcription 
factors belonging to 30 different families were affected 
under the drought condition. The genes encoding 
WRKY (2 genes; At2g40750; At2g23320), C2C2-CO-like 
(4 genes; At5g24930; At3g50410; Atlg30970; and 
Atlg43860), ARID (1 gene; Atlg76510), MADS (2 gene; 
At3g54340; At5g48670), EIL (1 gene; At3g20770), BZR (1 
gene; Atlg75080), and RAV (1 gene; Atlg68840) TFs were 
exclusively affected in sensitive genotypes (RAHS-14 and/ 
or RAHS-IPS-187 (Figure 7). The most attention gainers 
were MADS, ARID, and CO-like TFs, whose related genes 
are associated with plant development and periodic flow- 
ering. Other TFs whose encoding genes were exclusively 
present in Vagad/or GujCot-21 include ZF-HD (1 gene; 
Atlg74660), CPP (1 gene; At4gl4770), TCP (1 gene; 
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Table 2 Expression level changes of various antioxidant genes under drought stress in cotton root 



*Contigs represents antioxidant genes 


*TPM Count/tFold change 


Discription 






GUjCot-21 


RAHS-IPS187 




contig01875 




6.753 


- 


Cu/Zn superoxide dismutase 


contig00529 




2.996 


- 


superoxide dismutase-like protein [Oryza sativa Japonica Group] 


contig00696 




2.942 


- 


superoxide dismutase 


contig01057 




2.826 


- 


Superoxide dismutase [Mn]; Flags: Precursor 


contig01529 




2.454 


- 


copper chaperone for superoxide dismutase 
[Brachypodium distachyon] 


contig01075 




2.41 1 




Superoxide dismutase [Mn]; Flags: Precursor 


contigOl 149 




2.41 1 




Cu/Zn superoxide dismutase 


contig00854 




2.826 




, i • i , -I i-/-— • i • , -i 

cytosohc ascorbate peroxidase [Gossypium hirsutum] 


contig00876 




2.106 




cytosolic ascorbate peroxidase 1 [Gossypium hirsutum] 


contig00801 




2.904 


- 


cytosolic ascorbate peroxidase 1 [Gossypium hirsutum] 


contigOl 71 6 




3.148 


- 


Glutathione S-transferase 8 


contig02159 




2.826 


- 


Glutathione S-transferase 8 


contigOl 551 




2.826 


- 


Glutathione S-transferase 8 


contig02357 




2.826 




glutathione peroxidase [Populus trichocarpa] 


contigOl 51 8 




2.411 


- 


Glutathione S-transferase 8 


contig00449 




2.367 




Glutathione S-transferase 8 


contig00523 




2.205 




Glutathione S-transferase 8 


contigOl 789 




2.826 


- 


cytochrome P450 H202-dependent urate-degrading 
peroxidase [Glycine max] 


contig021 91 




2.826 




peroxidase [Ectocarpus siliculosus] 


contig00287 




2.744 


- 


PEROXIDASE 1 1 


contig00850 




2.563 


- 


cytochrome P450 H202-dependent urate-degrading 
peroxidase [Glycine max] 


contig00752 




2.563 




cytochrome P450 H202-dependent urate-degrading 
peroxidase [Glycine max] 


contig00393 




2.411 




putative peroxidase [Oryza sativa Japonica Group] 


contigOl 220 




2.411 




cytochrome P450 H202-dependent urate-degrading 
peroxidase [Glycine max] 


contigOl 737 




2.411 


- 


PEROXIDASE 1 1 


contigOl 277 




2.241 


- 


cytochrome P450 H202-dependent urate-degrading 
peroxidase [Glycine max] 


tProbe set represents antioxidant genes 


Vagad 


RAHS-14 




GhiAffx.53252.1.S1_ 


_s_at 




5.7884912 


glutamate synthase 


GhiAffx.1589.38.A1. 


_at 




4.280462 


SGT1 A 


GhiAffx.3996.1.S1_s 


_at 




3.7647517 


glutamate decarboxylase 


GhiAffx.4043.6.S1_at 




3.3622708 


glutamate decarboxylase 5 


GhiAffx.40354.1.S1_ 


.at 




3. 1 3/ buoy 


glutamate decarboxylase 5 


GhiAffx.30835.1.S1_ 


_s_at 




2.156465 


Peroxidase 1 


GhiAffx.221 15.1.S1_ 


x_at 




2.1306682 


Glutathione S-transferase 


GhiAffx.36353.1.S1_ 


.at 


5.4907985 




glutathione transferase 


GhiAffx.62078.1.S1_ 


.at 


3.2424738 




Glutaredoxin C9 


GhiAffx.21 573.1. S1_ 


.at 


2.9193442 




Peroxidase 10 


Ghi.6564.1.S1_at 




2.6120005 




peroxidase 


Ghi.6559.1.S1_at 




2.31927 




Peroxidase 73 


Ghi.6746.1.A1_s_at 




2.008061 




peroxidase 
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Table 3 Differentially expressed genes involved in membrane transport channels and cell wall synthesis under drought 
stress in Vagad and RAHS-14 



(a) Membrane transporter genes family 


Genotypes 


Cotton probe set ID 


TAIR 10 hit 


Fold change 


p-value 


Descriptions 


RAHS-14 0 




Ghi.1 932.1 .S1_at 


AT3G46740 


6.862507 


0.0339456 


outer membrane translocase 




Ghi.3760.1.S1_s_at 


AT1G19910 


2.9521239 


0.02668785 


ATPase, F0/V0 complex, subunit C protein 




Ghi.4450.1.A1_at 


AT3G01390 


2.2963412 


0.01776607 


V-type proton ATPase subunit G1 




Ghi.5565.6.S1_s_at 


AT1G15520 


2.354963 


0.02798233 


ABC transporter G family member 40 




Ghi.9771.1.S1_x_at 


AT2G33040 


3.4228158 


0.02484926 


Mitochondorial H + - or Na + -translocating F-type 




GhiAffx.5294.2.S1_x_at 


AT2G42210 


3.0888186 


0.03420474 


inner membrane translocas 




Gra.3081.1.A1_s_at 


AT4G32530 


8.158294 


0.04811877 


V-type and A-type ATPase (F-ATPase) Superfamily 




Gra.700.1.A1_s_at 


AT1G 15690 


3.1466656 


0.00362699 


vacuolar membrane proton pump 


RAHS-14 C 




Ghi.3763.1.A1_s_at 


AT3G63380 


5.7281384 


0.00042 


The P-type ATPase (P-ATPase) Superfamily 


Vagad 0 




Ghi.8876.1.A1_at 


AT1G51650 


2.1752877 


0.0222851 


Mitochondorial H + - or Na + -translocating F-type 




Ghi.9195.1.A1_s_at 


AT1G78720 


3.98914 


0.0432709 


preprotein translocase Sec61 alpha subunit 




GhiAffx.1 3955.1. A1_at 


AT3G47730 


2.8928177 


0.02658764 


The ATP-binding Cassette (ABC) Superfamily 




GraAffx.1 0962.1. A1_at 


AT3G47780 


2.3352501 


0.04103383 


ABC transporter A family member 7 




GraAffx.1 8588.1. A1_at 


AT3G59140 


2.7955472 


0.04233996 


The ATP-binding Cassette (ABC) Superfamily 


Vagad c 




GraAffx.29143.1.S1_s_at 


AT3G 13080 


2.0128014 


0.0003014 


The ATP-binding Cassette (ABC) Superfamily 


(b) Cell wall synthesis related gene family 


RAHS-14 D 




Ghi.3798.1.A1_s_at 


AT2G21770 


5.2988324 


0.02936104 


cellulose synthase 




Ghi.6435.1.A1_at 


AT5G24090 


12.502564 


0.02899352 


acidic endochitinase (CHIB1) 




Ghi.653.1.S1_s_at 


AT5G05170 


2.3362005 


0.04479192 


cellulose synthase (Ath-B) 




Ghi.825.2.S1_x_at 


AT2G40610 


7.838081 


0.00000635 


expansin, putative (EXP8) 




GhiAffx.1 589.84.S1_x_at 


AT5G17310 


2.6701603 


0.02671632 


UTP-glucose-1 -phosphate uridylyltransferase 




GhiAffx.1 81 56.1 .S1_at 


AT3G55250 


2.8398817 


0.04666388 


expressed protein 




GhiAffx.3601 8.1 .51 _at 


AT4G23920 


2.1763752 


0.01826835 


UDP-glucose 4-epimerase 




GhiAffx.46297.1.S1_s_at 


AT4G 17030 


3.4000947 


0.02354693 


UTP-glucose-1 -phosphate uridylyltransferase 




GhiAffx.49705.1.A1_s_at 


AT1G26810 


3.7181776 


0.048232 


galactosyltransferase family protein 




GhiAffx.52955.1.A1_at 


AT4G26940 


2.0691934 


0.0303665 


galactosyltransferase family protein 




GhiAffx.5996.1 .S1_s_at 


AT3G29090 


3.9073741 


0.02458842 


galactosyltransferase family protein 




GraAffx.29883.2.S1_s_at 


AT4G39210 


10.707062 


0.00723294 


ADP-glucose pyrophosphorylase 


RAHS-14 C 




Ghi.1 0493.1. S1_s_at 


AT5G57560 


8.634004 


0.0000114 


xyloglucan:xyloglucosyl transferase (TCH4) 




Ghi.1 127.2.A1_x_at 


AT4G25810 


7.3133554 


0.00000863 


xyloglucan:xyloglucosyl transferase (XTR6) 




Ghi.1 272.1. A1_s_at 


AT2G28470 


2.0890183 


0.0000454 


beta-galactosidase 




Ghi.4531.1.A1_at 


AT4G20460 


2.1103513 


0.00137449 


NAD-dependent epimerase 




GhiAffx.7421 .1 .S1_s_at 


AT5G44030 


50.340942 


0.0000394 


cellulose synthase (IRX5) 




Gra.1 335.1. A1_s_at 


AT5G22740 


2.0989435 


0.00195653 


glycosyl transferase family 


Vagad D 




Ghi.1 127.4.S1_x_at 


AT4G25810 


2.0596108 


0.02497086 


xyloglucan:xyloglucosyl transferase 




Ghi.5287.1.S1_s_at 


AT4G34150 


3.1211743 


0.03486373 


pectinesterase family protein 




Ghi.6619.1.S1_at 


AT2G39700 


2.451212 


0.04148584 


expansin, putative (EXP4) 




Ghi.6749.1.S1_s_at 


AT4G38400 


2.9026875 


0.04047799 


expansin family protein (EXPL2) 
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Table 3 Differentially expressed genes involved in membrane transport channels and cell wall synthesis under drought 
stress in Vagad and RAHS-14 (Continued) 
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cApdl Ibll 1, pULdLIVc ^LArDj 


Gra.1 984.1. A1_at 


AT1G04680 


4.803666 


0.00806574 


glycoside hydrolase family 


Gra.2243.1.A1_s_at 


AT5G09760 


2.1072824 


0.04628986 


pectate lyase family protein 


Gra.2265.1.A1_at 


AT5G 13870 


2.4477193 


0.04345882 


xyloglucan:xyloglucosyl transferase (EXGT-A4) 


Gra.2300.1.A1_at 


AT2G45470 


4.5964513 


0.03349594 


pectinesterase family protein 


Gra.2380.1.A1_s_at 


AT1G67750 


2.7641332 


0.01198629 


glycoside hydrolase family 


Gra.2812.1.S1_s_at 


AT4G03210 


2.5062635 


0.02250669 


xyloglucan:xyloglucosyl transferase 


Gra.377.1.S1_s_at 


AT2G06850 


3.0282483 


0.02482058 


xyloglucan:xyloglucosyl transferase (EXGT-A1) 


GraAffx.8742.1.S1_at 


AT1G32170 


15.896784 


0.00978865 


xyloglucan:xyloglucosyl transferase (XTR4) 


Vagad c 


Ghi.5772.1.S1_a_at 


AT4G 17030 


5.8524 


6.0179E-07 


expansin-related 


Ghi.6435.1.A1_at 


AT5G24090 


2.62142 


1.7468E-06 


acidic endochitinase (CHIB1) 


GraAffx.8742.1.S1_at 


AT1G32170 


3.4125 


0.00012776 


xyloglucan:xyloglucosyl transferase (XTR4) 


Ghi.1 0826.1. S1_at 


AT3G57270 


28.2959 


2.4257E-07 


glycosyl hydrolase family 



RAHS-14 D and Vagad 0 shows probe sets 
Arabidopis homologous ID at criteria 50 c , 



analysis under drought stress and RAHS-14 C and Vagad c under control condition. Cotton Probe converted unto 
x o query match and e-value cutoff 1e- 10 . 



At3g27010), REM (1 gene; At4g33280), G2-like (1 gene; 
At4gl3640), CAMTA (1 gene; At5g64220), and ABB VP 1 
(1 gene; At3gl8990). The ABB VP 1 has been designated 
as being heat and drought responsive. The CAMTA func- 
tions as a negative regulator under biotic stress and is 
described as being a cold responsive gene. Similarly, the 
induction of G2-like TF is said to confer thermo-tolerance 
to the plants. The bZIP (6 genes; Atlg77920; At3gl9290; 
At4g34590; At2g40950; At2g42380; and At5g28770) and 
ARF (3 genes; At4g23980; At5g60450; and At5g62000) 
were predominant in the Vagad/or GujCot-21, which are 
reported to be regulated by phytohormones [39,40] .There 
were few genes encoding various TFs whose expression 
was almost unanimous in all four genotype, such as MYB, 
NAC, and AP2/EREBP. These TFs are well reported as 
playing roles in developmental stages. The MYB TF is 
required for fiber development [41]; the NAC is a plant- 



specific transcription factor that is known to play a varied 
role in plant developmental processes [42]. 

Tissue-specific expression of genes in cotton root may 
contribute to drought tolerance 

In order to identify the tissue-specific genes in roots that 
respond to drought stress, differentially expressed genes 
in transcriptome analysis were carefully curated for their 
corresponding homologues in Arabidopsis, and these 
Arabidopsis IDs were queried against Genevestigator 
(Figure 8 A and B) [43]. The differentially expressed 
genes showed specific expression in different root zones, 
including root hair, lateral root, elongation zone, meris- 
tamatic zone, stele, and pericycle as well as in the root 
tip (Figure 9, Figure 8 A and B). Genevestigator results 
showed that the expressions of genes, such as ERD 
(early-responsive to dehydration stress), family protein 
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Figure 7 Differentially expressed TFs in all genotypes (Vagad, GujCot-21, RAHS-IPS-187 and RAHS-14). 
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(AT4G22120), ABRE binding factor (AT3G19290), tono- 
plast intrinsic protein (AT2G36830), ribosomal protein 
(AT1G17080), 26S proteasome (AT4G24820), methyl- 
transferase (AT4G34050), ADH, alcohol dehydrogenase 
(AT5G43940), and so on, were higher in root hair of the 
tolerant genotype (Figure 9, Figure 8A and B); whereas in 
RAHS-14 genes, beta-fructofuranosidase (AT1G35580), 
kinases-related proteins (AT3G48530), cellulose synthase 
(AT4G39350), and unknown protein (AT5G02020) 
showed higher expression. Similarly, the lateral root zone 
showed the expression of only one gene that exhibited 
oxidoreductase activity (AT1G60710) in tolerant genotype 
and two genes encoding vacuolar membrane-related pro- 
tein (AT3G26810) and signal transduction histidine kinase 
protein (AT3G21510) in sensitive genotype (Figure 9). In 
tolerant genotype, a greater number of genes were 
expressed in the elongation zone of the root as compared 
with sensitive genotype. The genes encoding tonoplast 
intisic protein (AT1G09330), cytosolic phosphoglucomu- 
tase (AT1G23190), ketoacyl-CoA synthase (AT2G26640), 
26S proteasome (AT4G24820), and cellulose synthase 
(AT5G05170) were expressed in tolerant genotype. Simi- 
larly, the genes encoding phosphogluconate dehydrogen- 
ase (AT1G64190), aldehyde dehydrogenase (AT1G74920), 
and tetratricopeptide repeat (TPR)-like superfamily pro- 
tein (AT4G 10840) were expressed in sensitive genotype 
(Figure 9, Figure 8A and B). Interestingly, tolerant geno- 
type expressed a large number of genes in meristematic, 
pericycle, and stele regions of the root as compared 
with sensitive genotype. Stele and pericycle regions of 
the tolerant root expressed genes such as transcrip- 
tion factor, FARl-related (AT1G10240), transmem- 
brane nitrate transporter (AT1G32450), coumarate: CoA 
ligase (AT1G51680), Aldo/keto reductase (AT1G60710), 
receptor-like protein containing leucine-rich repeats 
(AT1G65380), zinc finger protein (AT1G74660), and so 
on; whereas the genes expressed in the pericycle of sensi- 
tive genotype were related to the ethylene-responsive 
family protein (AT4G29100), DAHP (3-deoxy-7-phospho- 
heptulonate synthase) synthetase family protein, related to 
senescence (AT1G22410), amino-acid transporter family 
protein (AT1G47670), and alcohol dehydrogenase 
(AT1G74920) (Figure 9, Figure 8A and B). The expression 
analysis of the differentially expressed genes by Genevesti- 
gatotr revealed that approximately 44 and 10 genes were 
expressed in the meristematic tissue of tolerant and sensi- 
tive genotype, respectively. In tolerant genotype, many 
genes related to ribosomal family protein (AT3G06680, 
AT2G34480, AT3G16780, AT5G22440, AT1G27400, 
AT2G19740, AT3G53740, AT3G09630, AT5G02870, 
AT4G10450, AT3G62870, AT3G47370, AT4G33865, and 
AT4G34555), NAPl-related protein 1(AT1G74560), chro- 
matin modifying enzymes (AT5G64610), and importin pro- 
teins (AT3G06720) were present. Similarly, in sensitive 



genotype genes related to alcohol dehydrogenase 
(AT1G74920), ribosomal protein (AT2G44065), calmodulin 
binding protein (AT5G20720), and protein kinases 
(AT5G47750) were highly expressed (Figure 8B). The pres- 
ence of the NAPl-related transcript at a meristematic level 
serves as a good indicator of our root zone profiling of the 
differentially expressed genes. At the root tip region of the 
tolerant genotype, we observed the expression of the genes 
related to ATP binding and aminoacyl-tRNA ligases 
(AT1G09620), ribosomal family protein, (AT2G34480, 
AT3G09630, AT3G47370, and AT3G53740), and cytosolic 
glucose-6-phosphate dehydrogenase (AT5G40760); whereas 
the expression of these differential genes was absent in sen- 
sitive genotype (Figure 8A). 

Discussion 

The present study reports the genome-wide expression 
studies conducted on the roots of drought-tolerant and 
sensitive genotypes of G. herbaceum and identifies the 
various genes and pathways responsible for drought tol- 
erance. The drought- tolerant (Vagad and GujCot-21) 
and -sensitive genotypes (RAHS-14 and RAHS-IPS-187) 
showed significant variations in root structure and their 
length under control and mannitol stress (Figure. 1). 
The longer root length in Vagad and GujCot-21 geno- 
types showed their inherent ability to adapt to the water 
stress condition [44]. The adaptive advantage of 
increased root length and more number of root hairs 
that facilitate access to water under drought stress has 
been reported earlier [45]. Vagad and GujCot-21 dis- 
played hydrotrophism at variable concentrations of man- 
nitol, in which the root detects a stress gradient and 
directs its growth [46]. In Gujcot21 and Vagad, longer 
primary roots at 6% and 8% mannitol reflected an inter- 
esting plant adaptation for better survivability (Figure 1). 
At moderate stress of 4% mannitol, Vagad and GujCot- 
21 had improved anatomical adaptation by smaller leaf 
size and inward rolling of the leaves in order to restrict 
the exposed leaf surface, thereby limiting water loss by 
evaporation (Figure 1) [47]. The expression profiling by 
microarray analysis and transcriptome sequencing of the 
root revealed that a large number of genes such as 
osmotin, serine, aldo-keto reductase, LEA, chalcone syn- 
thase, chitinase, RD29, RD19, proline oxidase, ERD9, 
sodium hydrogen exchanger, glycin, alcohol dehydrogen- 
ase, brassinosteroids, HSP 70, and metalothionin were 
up-regulated in Vagad and GujCot-21 (Additional file 4, 
Additional file 5, and Additional file 6). Many of these 
genes were reported to help in osmoregulation and, 
thus, protect plants in stress management [4,10,29,38,48]. 
Many of the differentially expressed genes identified 
were found to represent diverse metabolic pathways 
(Additional file 4, Additional file 5 and Additional file 6). 
It was interesting to note that sensitive and tolerant 
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Figure 9 Differentially expressed genes (probe sets and contigs) analyzed by Genevestigator in mapping the specific expression of 
genes in different root zones. 



genotypes underwent a larger degree of transcriptional 
modulation representing the various metabolic processes 
involved during drought stress (Figure 2, Figure 3 and 
Figure 5). This transcriptome modulation may eventually 
result in the synthesis of the secondary metabolic com- 
pounds that provide adaptive advantages to genotypes 
[48]. The gene ontology analysis of the differentially 
expressed genes obtained in microarray and transcrip- 
tome sequencing expression indicates a cross talk 
among various pathways that helps the tolerant geno- 
types to survive under drought stress. Many GO terms 
were exclusively present in Vagad/or Gujcot-21, and 
most of them were related to categories of abiotic 
stress, such as cellular response to stimulus, response 
to hormone stimulus, response to abiotic stimulus, ion 
transporter, kinase activity (Figure 3), response to inor- 
ganic substances, response to water deprivation, re- 
sponse to jasmonic acid, response to salinity, response 
to hyperosmotic stress, and response to salinity 
(Figure 5). The enrichment of these stress-related GO 
terms are indicative of better management of drought 



stress in tolerant genotypes [26]. The GO terms such 
as reactive oxygen species, lipid localization, lignin, su- 
crose, glucan, and starch biosynthesis processes were 
enriched in RAHS-14/or RAHS-IPS-187, but the per- 
centages of stress-inducible GO terms and genes were 
less as compared with Vagad/or GujCot-21 (Figure 2 
and 5). The result probably indicates that sensitive gen- 
otypes had a slow response and, hence, gradually pre- 
pared themselves to adapt to drought stress; whereas 
drought-tolerant genotypes were already equipped with 
these processes, therefore had a prompt response to- 
ward stress, and, hence, adapted better to stress. The 
KOBAS analysis showed significant enrichment of beta- 
alanine and histidine metabolism and flavanoid metabol- 
ism in GujCot-21 (Figure 6A), which indicates the 
formation of osmoprotective molecules under the 
drought condition. During drought stress, beta- alanine 
and histidine have been shown to increase in several 
different plant species, including Arabidopsis, maize, and 
wheat [49-51], and are considered as having an osmopro- 
tective function. The secondary metabolites, such as lignin 
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precursors, and flavonoid increase in plants under stress 
[52]. The presence of ABA phytohormonal pathways in 
RAHS-IPS-187 (Figure 6B) could be concurrent with 
their response to drought stress by mediating and 
channelizing many stress-responsive genes that help the 
plants in their survival over stress [53]. In contrast, 
Vagad and GujCot-21 adjust osmotic balance under 
drought stress by an increase in osmoticum deposition 
such as beta-alanine and catabolism, import, and 
utilization of various other secondary metabolites 
(Figure 6A, Additional file 5 and Additional file 6). The 
extent of osmotic adjustment in the root, although sub- 
stantial, is insufficient to maintain turgor wall pressure 
under drought stress [48,54,55]. The transporter pro- 
teins and ion channels also play an important role in 
osmoregulation or drought tolerance [50]. Various 
membrane transporters, including ion (ABC transpor- 
ters, ion channels) or sugar transporter, were enriched 
in Vagad under drought stress (Table 3). These trans- 
porters could both maintain the homeostasis of intracel- 
lular ion concentrations of cells and stabilize the 
physiological balance of plants under drought stress 
[29]. The higher expression of more number of genes 
related to membrane transporters under drought stress 
in tolerant genotypes (Vagad and GujCot-21) showed 
their plastic or adaptive behavior toward drought stress 
(Additional file 6, Table 3). Cell wall-loosening proteins 
are believed to play key roles in controlling cell wall ex- 
tension [56]. The expression of expansins and xyloglu- 
can endotransglycosylase (XET) genes were found to be 
higher in Vagad when drought was imposed (Additional 
file 5). It was reported earlier that expansin activity and 
extractable expansin protein increased in the root when 
plants were exposed to drought stress [57]. The higher 
expression of expansin under drought stress indicates 
changes in cell wall structure or the chemistry of toler- 
ant genotypes that probably facilitate turgor pressure 
and the maintenance of ionic balance [56]. Genes en- 
coding for peroxidases showed a significantly higher ex- 
pression of ascorbate peroxidase, SOD in Vagad and 
GujCot-21 (Table 2). The expression of these oxidative 
radical scavenging enzymes indicates the better manage- 
ment of oxidative radicals in tolerant genotypes during 
drought stress [34]. Drought resulted in elevated ex- 
pression levels of several genes in RAHS-14 and 
RAHS-IPS-187 like 1-aminocyclopropane-l-carboxylate 
oxidase (Ghi.714.1.Sl_s_at; in Additional file 3 and 
Additional file 6), which catalyze the final step of 
ethylene biosynthesis [58]. Ethylene is elicited and 
acts as a signal mediator under almost all biotic and 
abiotic stress conditions, including drought [58]. The 
expression of ERF-related genes in sensitive genotypes 
under the drought condition indicates that sensitive 
genotypes responds to drought mainly by the ethylene 



pathway, leading to senescence [59]. The results were 
in agreement with our previous results on leaf tran- 
scriptomes for these genotypes [31]. Among TFs, the 
expression of the transcription factors containing the 
AP2-EREBP domain were higher in tolerant genotypes 
under drought stress (Figure 7); similarly, bHLH (basic 
helix loop helix) and bZIP (basic leucine zipper domain) 
transcription factor in Vagad and GujCot-21 were 
higher. TFs bHLH and bZIP belongs to a family of pro- 
teins that are specific to plants and are found to play a 
role in a diverse set of biological processes under biotic 
and abiotic stress. In case of Oryza sativa, the over- 
expression of bHLH and bZIP transcription factors 
resulted in the expression of several stress-responsive 
genes in the transgenic plants, thereby conferring the 
drought tolerance [53,60]. The transcription factors 
WRKY 75 (GhiAffx.6177.1.Sl_at, GhiAffx.61 976.1. 
Sl_at) and ZAT 10 (contig00573) were exclusively 
expressed in Vagad/or Gujcot-21 and are known to 
regulate root hair development and play a role in 
stress conditions [17]. It was reported that RHD3 
(Root hair defective 3) and LBD (Lateral organ 
Boundaries Domain) genes play important roles in 
root development and the formation of root hairs in 
Arabidopsis [10,60,61]. The higher expression of 
RHD3 and LBD genes in tolerant Vagad and GujCot- 
21 could be considered the key regulatory factors for 
better root length development and for maintaining 
the plasticity of roots under drought stress conditions 
(e.g. osmotic adjustments in roots) (Additional file 5 
and Additional file 6). There was a preferential expres- 
sion of genes, such as ERD (early responsive to dehydra- 
tion stress) family protein, tonoplast intrinsic protein, 
ribosomal protein, and methyltransferase in the root hair; 
oxidoreductases in the lateral root region; and cellu- 
lose synthase, tonoplast intisic protein, and cytosolic 
phosphoglucomutase in the elongation region of the 
root in Vagad (Figure 8 and 8A). Cell division-related 
genes encoding NAP1 proteins, ribosomal proteins, and 
transcription factors were over-represented in the meri- 
stem. It should be noted that preferentially expressed 
genes in the pericycle and stele regions of the Vagad 
were related to nitrate transporter, coumarate-CoA lig- 
ase, and aldo/keto reductases (Figure 9 and 8A). Nitrate 
is the major nitrogen source for the synthesis of amino 
acids and nucleic acids in plants. In addition to being 
assimilated in the cytoplasm, it has been proposed that 
nitrate may be an important osmotic solute [62]. The 
preferential expression of the nitrate transporter in 
Vagad may help in the cellular homeostatic mainten- 
ance of plants under drought stress. We also observed 
that many genes which encode cell wall modification 
enzymes, such as cellulose synthase and cell expansin 
in the lateral root region of Vagad, are considered as 
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Description: Excel file containing the list of unique up regulated genes 
and their annotation of RAHS-14 under drought stress. 

Additional file 4: Annotation and fold change (fold change >2) of 
up regulated genes (probe sets) in Vagad under control condition. 

Description: Excel file containing the list of unique up regulated genes 
and their annotation of Vagad under control condition. 

Additional file 5: Annotation and fold change (fold change > 2) of 
up regulated unique genes in Vagad under drought stress. 

Description: Excel file containing the list of unique up regulated genes 
and their annotation of Vagad under drought stress. 

Additional file 6: Contigs analysis by R-statistics of GujCot-21 and 
RAHS-IPS-187 genotypes. Description: Excel file containing detail 
analysis of contigs in both genotypes. Worksheet 1 containing total 
differential contigs, worksheet 2 is up-regulated contigs in RAHS-IPS-187, 
worksheet 3 is up-regulated contigs in GujCot-21, worksheet 4and 5 are 
agriGO of RAHS-IPS-187 and GUjCot-21 respectively. Worksheet 6 and 7 
are KOBAS analysis of RAHS-IPS-187 and GUjCot-21 respectively. 

Additional file 7: BLASTX analysis of contigs and singlets against 
NCBI NR database. Description: Excel file containing BLASTX (TAIR 
database) results of supercontigs (worksheetl), super singleton 
(worksheet2), contigs of GujCot-21 (worksheet3), singletone of GujCot-21 
(worksheet^, contigs of RAHS-IPS-187 (worksheet 5) and singletone of 
RAHS-IPS-187 (worksheet 6). 

Additional file 8: BLASTX analysis of both transcriptome contigs 
and singlets against TAIR database. Description: Excel file containing 
BLASTX (NCBI NR database) results of supercontigs (worksheetl), super 
singleton (worksheet2), contigs of GujCot-21 (worksheets), singletone of 
GujCot-21 (worksheets, contigs of RAHS-IPS-187 (worksheet 5) and 
singletone of RAHS-IPS-187 (worksheet 6). 

Additional file 9: BLASTN analysis of both transcriptome contigs 
and singlets against publicly available cotton EST sequences. 

Description: Excel file containing BLASTN (Cotton EST) results of 
supercontigs (worksheetl), super singleton (worksheet2), contigs of 
GujCot-21 (worksheet3), singletone of GujCot-21 (worksheet^, contigs of 
RAHS-IPS-187 (worksheet 5) and singletone of RAHS-IPS-187 
(worksheet 6). 

Additional file 10: Venn diagram of contigs obtained after BLASTX 
in NCBI-NR, TAIR and Cotton EST database (BLASTN). Contigs of (A) 
Drought tolerant genotype (GujCot-21) (B) Sensitive genotype 
(RAHS-IPS-187) (C) Supercontigs Description: PPT contains a figures of 
venn diagram representing unique and common contigs in different 
database. 



playing a key role in border cell separation and the 
elongation of the root [55-57]. Our results showed that 
a significant number of genes encoding ribosomal pro- 
teins were highly expressed (fold change >2) in tolerant 
genotypes (Figure 9 and Additional file 6). Plants 
expressed a high number of ribosomal proteins to bal- 
ance cellular protein synthesis in response to environ- 
mental variations and, at the same time, to adapt to the 
environment [63]. Earlier studies showed that genes en- 
coding ribosomal proteins play an important role in the 
development of root meristem and contribute to the 
growth and development of the root system [63-65]. A 
higher expression of ribosomal protein could be consid- 
ered important during cell growth and proliferation 
under drought stress in the case of tolerant genotypes. 
Differentially expressed genes were analyzed using the 
Genevestigator tool (https://www.genevestigator.com/gv/ 
plant.jsp) (Figure 8A and B) for their specific expression 
in different root zones [43]. Interestingly, a large num- 
ber of these genes belonged to ribosomal coding proteins 
and were expressed in either the meristematic zone or the 
root tip region, which could maintain the homeostatic bal- 
ance of protein synthesis of the plant under the stress con- 
dition [66]. This comprehensive analysis of the 
transcriptional profile during drought stress will advance 
our fundamental understanding of the various genes and 
major metabolic pathways that provide direction for the 
future genetic engineering of drought-stress tolerance in 
cotton. 

Conclusion 

The present study highlights the impact of drought on 
the stress-responsive signaling pathways in roots, 
which regulates the plant adaptation to limited water 
condition as well as other environmental stresses. The 
stress-induced transcripts had a variety of membrane 
transporters, such as peroxidase, catalase, WRKY, AP2 
transcription factor, stress-responsive genes, and meta- 
bolic and regulatory pathways. The root transcriptome 
analysis revealed several important genes regulating 
drought tolerance and that could be considered prob- 
able potential targets for candidate gene selection in 
improving the tolerance of plants. 

Additional files 



Additional file 1: List of primer used in RT-PCR Description: Word 
file containing list of primer sequences used for RT-PCR. 

Additional file 2: Annotation and fold change (fold change > 2) of 
up regulated unique genes in RAHS-14 under control condition. 

Description: Excel file containing the list of unique up regulated genes 
and their annotation of RAHS-14 under control condition. 

Additional file 3: Annotation and fold change (fold change > 2) of 
up regulated unique genes in RAHS-14 under drought stress. 
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